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ABSTRACT 

A  program  is  developed  to  ascertain  the  response  of  a 
frictionless  water  column  to  internal  waves.   Using  density 
strata  as  input,  the  program  selects  ten  frequencies  at 
equal  intervals  within  a  spectrum  of  internal  waves  bounded 
above  by  the  maximum  Vaisala  frequency  and  below  by  the 
inertial  frequency.   Ray  paths  within  the  column  are  plotted 
for  these  frequencies.   The  first  ten  normal  modes  for  each 
frequency  are  computed.   The  first  three  modes  are  plotted. 
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I.   INTRODUCTION 

Internal  waves  exist  in  all  oceans,  most  bays  and  lakes; 
they  vary  widely  in  amplitude,  period,  and  depth  (La  Fond, 
1962).   Recently,  advances  have  been  made  in  interpreting  the 
major  features  of  the  oceanic  internal  wave  field  (Vunsch  and 
Dahlen,  1970). 

Internal  waves  can  be  observed  as  density  perturbations 
through  a  stable,  stratified  water  column.   Since  temperature 
correlates  with  density  in  the  majority  of  the  upper  ocean 
areas,  internal  waves  often  appear  as  oscillations  about  the 
region  of  the  thermocline,  exhibiting  greatest  amplitude  and 
phase  speed  where  Ine  tnermocliRe  is  tiie  most  intense.   J-iie 
oscillation  initiated  within  the  water  column  is  described  in 
terms  of  the  Vaisala  or  stability  frequency.   This  frequency 
may  be  derived  as  follows  (Zalkan,  1966). 

Consider  the  force  acting  on  a  water  parcel  raised 
vertically  from  its  equilibrium  position: 

F  =  g[e-(e+r^Az)]dV-p^Azdv        (1) 

0  c 

where:   g = gravitational  acceleration, 

r  =  density , 
&Z  -  vertical  displacement, 
dV = minute  volume  of  the  parcel,  and 
c  = sound  velocity.   The  first  term  on  the  right  hand 
side  is  the  restoring  force  due  to  buoyancy.   The  second  term 
is  the  contributing  force  of  compressibility.   The  equation 


for  the  induced  vertical  motion  is: 

,2  -no 


edV^=-gSZdV-^ zdv  <2> 

The  harmonic  solution  to  equation  (2)  is  the  Vaisala  frequency 

defined  as 

/   aO    2 

(3) 

For  oceanographic  applications  the  water  is  effectively  in- 
compressible and  the  second  term  under  the  radical  is 
negligible.   Thus, 

N(Z)  =  /^S  (4) 

The  Vaisala  frequency  is  greatest  in  magnitude  where  the 
thermocline  is  most  intense  (see  Figure  l). 

The  vertical  motion  initiated  by  the  internal  wave  pac- 
ing through  the  water  column  is  governed  by  the  differential 
equation  (Fjeldstad,  1933): 

p<?i>+0,"T^(N2"tf2)¥=°  <5> 

where:   V = vertical  water  velocity, 

k = internal  wave  number, 

(f  =  internal  wave  frequency,  and 

f  = inertial  frequency.   Equation  (5)  is  a  hyperbolic 
wave  equation  whose  characteristics  describe  the  direction 
of  propagation  of  internal  wave  energy.   Vith  the  appropriate 
boundary  conditions  that  no  vertical  flow  exists  at  the  sur- 
face or  at  the  bottom,  the  internal  wave  equation  satisfies 
the  criteria  for  a  Sturmian  system.   The  imposed  boundary 
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VAISALA   FREQUENCY 


FIGURE  1.   Variation  of  Vaisala  Frequency 
with  Temperature  Profile 


conditions  are  satisfied  by  the  simultaneous  propagation  of 
two  waves  having  the  same  horizontal  wave  number,  one  with 
its  vertical  component  of  phase  velocity  up,  the  other  with 
its  vertical  component  down.   The  stipulation  that  these 
waves  cancel  at  the  boundaries  leads  to  the  structure  of 
normal  modes.   These  modes  constitute  a  discrete  set  of  wave 
numbers  that  are  the  eigenvalues  of  the  Sturmian  system. 

This  thesis,  using  as  a  basis  Fjeldstad's  classic  paper 
(Fjeldstad,  1933),  developed  a  numerical  program  that  assessed 
the  response  by  a  f rictionless ,  incompressible  water  column 
to  internal  waves  in  a  frequency  band  bounded  above  by  the 
maximum  Vaisala  frequency  and  below  by  the  inertial  frequency. 
Station  density  strata  was  used  as  input.   Frequencies  were 
selected  within  the  band  and  their  ray  paths  were  drawn. 
Also,  the  governing  differential  equation  was  numerically 
integrated  and  the  resulting  endpoints  were  refined  to  deter- 
mine the  first  ten  normal  modes.   The  first  three  modes  were 
plotted. 

The  program  was  first  tested  using  density  strata  with 
known  output.   The  program  then  treated  data  taken  from  ocean 
stations  off  the  coast  of  Monterey,  California  in  mid-August, 
1970. 
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II.   THEORY 

A.   DEVELOPMENT  OF  THE  INTERNAL  WAVE  EQUATION 

For  a  f rictionless ,  incompressible  water  column  the 
governing  equations  (following  Feldstad,  1933)  are: 
Momentum 

H  +  fU  +  lf|=0  (7) 


Continuity 


Densily  invariant 

&  +  "jl  +  *#  +  vg-o  <"» 

The  X-axis  is  directed  in  the  direction  of  wave  travel, 
the  Y-axis  is  in  the  horizontal  transverse  plane  of  travel, 
and  the  Z-axis  is  vertical,  positive  upward.   U,  V,  and  V  are 
the  velocity  components  of  the  water  in  the  X,  Y,  and  Z  axes, 
respectively.   f  is  the  coriolis  parameter.   P  is  the  pressure 
The  depth  is  assumed  to  be  uniform.   Z  =  0  is  the  bottom;  Z  =  h 
is  the  undisturbed  upper  surface.   For  a  steady  state 
condition , 


dU   dV   dV 
Thus, 


dt-dt-  Lt  =  0-  (11) 


dP  =  -r>gdZ.  (12) 
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P  can  be  represented  as 

p  =  po  +  g  J  ?dz'  (13) 

where  Pq  is  the  undisturbed  pressure  at  the  surface.   Let 

£=£0(Z)  +  ^(X.Y.Z.t),  (14) 

where  \~    is  the  undisturbed  density  and  ^,  the  perturbation. 
Therefore,  pressure  for  the  perturbation  is: 

P  =  P0+gJ  ^0dZ+Pl  (15) 

where 

(16) 


Jh 


Envoking  a  general  assumption  that  velocities  and  perturbations 
of  density  and  pressure  are  small  enough  that  their  products 
and  squares  may  be  neglected,  equation  (10)  may  be  written: 

If,  for  the  internal  wave  form,  harmonic  oscillations 

such  as  the  following  are  considered: 

u  ,v  , P  rv/  cos  (($  t  -  kx)  exp  (ay  ) 

V,V  r^  sin  (cf  t  -  kx  )  exp  (ay) 

where  the  constant  a  controls  the  transverse  envelope  and 

allows  for  the  possibility  of  edge  waves.   The  sign  of  a  is 

such  that  the  wave  amplitude  decays  away  from  the  boundary. 

The  governing  equations  now  become: 

Pn 

rr  u-  fv  +  k~  =  o  (18) 

p, 

(fV+fU+ag^-0  (19) 

^0 
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+  g+pL 

\  0 

J  +  aV  + 

dZ 

(20) 
(21) 


Solving  equations  (18)  and  (19)  for  U  and  V: 


nV*^!j  (22) 

V#^S-  (23) 


Solving  for  the  vertical  pressure  gradient: 


i  _n2  -qf2  q 


,V  (24) 


dz      cr'       ° 

2   -   ^  o 
where  N  =f)       An     as  previously  defined.   When  one  inserts 

(22)  and  (23)  in  (9) J 

£l  -  .  i£2+lh  M  (25) 

<^~  tf(k2-a2)  dZ  '  " 

Vhen  P,  is  eliminated  from  (24)  and  (25)  and  the  transverse 

boundaries  are  disregarded  (i.e.,  set  a  =  0), 

d2¥  ,    .  1     dH)xdV^  r_k2__irNT2     rf2-wr      n  /orx 

^2+(^lT)dZ+[072-772]tN     -0     ]V  =  0  (26) 

Equation  (26)  is  the  internal  wave  equation  describing  the 
response  by  a  water  column  to  a  plane  internal  wave  of  known 
frequency, v.   The  order  of  magnitude  of  the  second  term  is 
small  compared  to  the  first  and  third.   However,  the  numerical 
program  can  include  the  second  term  without  any  difficulty. 
It  will  be  so  retained. 
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B.   APPLICATION  OF  RAT  THEORY 

An  analytical  solution  to  equation  (26)  may  be  made  by 
the  transformation: 

V=  (Z)  exp  (^"  yf  §)•  (27) 

When  the  effect  of  the  second  term  in  eq.  (26)  is  small,  the 
analytical  solution  to  the  internal  wave  equation  can,  with 
a  high  degree  of  accuracy,  be  represented  by T  ,  the 
oscillatory  portion.   The  oscillatory  solution  is: 


Y 


+  (   k   o)(N2-fl2)W=0  (28) 


<r-f 

Solutions  to  equation  (28)  are  of  the  form: 

I  =A  exp  i  (k1X+k1  ^  2~     2   Y  ^^  ^29^ 

where  k-.  =  horizontal  wave  number. 

Define  0  as  the  angle  between  the  horizontal  and  the 
phase  velocity.   Then, 

k!                                tf2  -  f 2    1//2 
cos0  = ~ ~ T-p:  =  [— 5 x]  (30) 

(^     2.,  2rN2-^2lNi/2    N2  -  f 2 

(kl     +kl    [o^T^2J) 


Since 


where 


0        0    1/2  ?         ?    1/2  k 

(tf2-f2)  =  (N2-tf2)  2    1/2  (31) 

(k^+k/)172 

k   2N2  +  k   2f2    ^2 
<T     =[-3-5 1 ]  Ola) 

kl     +k2 

XT2         2    1/2 
^^V^  (32) 

|k|  =(k12  +  k22)l/2  (33) 
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Figure  2  illustrates  the  angle  between  the  horizontal 
and  the  phase  velocity.   The  phase  velocity,  broken  into  X 
and  Z  components,  is: 


X   ,  2  ,  ,  2~/,  2  _,_  ,  2x3/2  L  1 
12    ^  1     2  ' 


[k  2N2  +  k  2f2]l/2     (34) 


and 


C7= 2   2  2  Sh    [k  2N2  +  k  2f2]l/2 

Z   (k12  +  k22)3/2    1       2 


For  the  group  velocity, 


c„  =Vk  (T 


g 


where 


->     o        k  2N2  +  k  2f2    1//2 
n      -  V    I"  C_i -N 


g  Kl 


kl2  +  V 


(35) 


(36) 


(37) 


/'  o  o  \ 


The  X-component  of  group  velocity  is: 


(C    ) 


g'x~d 


A        k  2N2  +  k02f2    1/2 

[(-1— 5 V-)       J 


k 


2  2 

V  +  k2 


klV 


[■ 


k|3    (kx2N2 +  k2^f^) 


(N2-f2)       t    i 
2^2X1/2J 


The  Z-component  of  group  velocity  is: 


k  2N2  +  k2f2  ^2 

Z  0k2     k2  +  k2 


(39) 


kl2k2  L 

|k|3    (k12N2  +  k22f2)1 


(N2  -  f2) 


Til) 


(40) 


Define©   as  the  angle  made  with  the  horizontal  by  the 
group  velocity  vector.   Then, 
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FIGURE  2.   Illustration  of  Angle  between 
Horizontal  and.  Phase  Velocity 
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,    e     (Cg)z   "Ki    rN2-f2i1/2 

The  ray  trace  part  of  the  program  drew  the  ray  paths  as 
depicted  in  equation  (41)  for  ten  frequencies  selected  within 
the  inertio-gravity  band. 

Determining  the  dot  product  of  the  phase  and  group 
velocities : 


C  0     C 

g 

(kl(kl2N2  +  k22f2)l/2)(klkl2       (N2-f2)     } 
|k|3  |k|3  (k12N2  +  k22f2)1/2 

k0(k  2N2  +  k  2f2)l/2     k  2k0(N2-f2) 
_  (_£ i i )    ( 1 2 s 

|k|3  |k|3(k12N2  +  k22f2)1/2 

=  0  [A2) 

Thus  the  phase  velocity  is  perpendicular  to  the  group 
velocity.   Also, 

k2(k22-k12)f2 
Cz+(C^Z=(k12N2  +  k22f2)1/2(k12  +  k22)3/2        ^^ 

The  phase  velocity  and  group  velocity  then  have  a  relationship 
as  seen  in  Figure  3. 


C.   NORMAL  MODE  APPLICATION 

With  the  use  of  suitable  boundary  conditions,  equation 
(26)  may  be  solved  to  obtain  modal  structures.   The  applicable 
conditions  are: 

V=0    at    Z  =  0  (44) 
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z 


Figure  3.   Relation  between  Phase  and  Group  Velocities 
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and 

d!  ~  ^772  SW  =  °   at   Z  =  h  (44a) 

To  a  high  degree  of  approximation  the  vertical  velocity  shear 
is  negligible  at  the  surface  and  the  second  condition  becomes: 

V  =  0    at    Z  =  h  (45) 

Equations  (26),  (44),  and  (45)  compose  a  Sturmian  system. 
¥,  the  vertical  velocity,  is  the  eigenf unction.   For  a  given 
frequency,  discrete  wave  numbers  or  eigenvalues  are  determined, 
The  program  utilized  a  fourth-order  Runge-Kutta  numerical 
method  to  obtain  an  end  value  after  integration  of  the 
internal  wave  equation  over  the  water  column.   The  endpoint 
was  refined  using  a  Newton-Raphson  iteration  procedure. 
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III.   DEVELOPMENT  OF  THE  MODEL 

A.  MAIN  PROGRAM 

The  main  program  prints  the  station  density  stratifica- 
tion.  The  program  then  calls  upon  the  three  sub-programs  to 
interpolate  density  between  station  depths,  to  plot  ray  paths 
for  selected  internal  wave  frequencies,  and  to  print  out  and 
plot  modal  structures  for  these  frequencies. 

B.  SUBROUTINE  INTERP 

Subroutine  INTERP  interpolates  the  density  stratification 
for  400  incremental  depths  through  the  water  column.   The 
interpolation  is  accomplished  by  the  use  of  a  third-degree 
polynomial  that  matches  points  and  slopes  at  each  depth  where 
data  was  obtained  from  a  Nansen  bottle.   To  guarantee  density 
stability,  the  smaller  of  the  slope  above  or  below  each  point 
is  chosen  as  the  slope  at  that  point.   From  the  interpolation 
polynomial,  INTERP  computes  the  Vaisala  frequency  for  each 
density  value. 

C.  SUBROUTINE  RAY 

Subroutine  RAY  plots  the  interpolated  density  and  the 
Vaisala  frequency  as  a  function  of  depth.   RAY  then  determines 
the  inertio-gravity  frequencies  at  rates  of  l/lOO,  l/lO,  2/l0, 
3/10,  4/10,  5/10,  6/10,  7/10,  8/10,  and  9/10  of  the  frequency 
band.   The  program  then  plots  the  characteristic  paths  over 
the  depth  of  the  column  for  each  of  these  frequencies. 
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Additionally,  subroutine  RAY  determines  the  percentage  of 
the  water  column  depth  where  the  selected  frequency  exceeds 
the  Vaisala  frequency.   This  percentage  represents  that  por- 
tion of  the  column  where  the  Sturmian  solution  to  the  internal 
wave  equation  is  exponential.   A  greater  percentage  implies 
increasing  difficulty  in  applying  boundary  conditions  during 
the  numerical  integration.   The  percentage  figure  for  each 
frequency  is  listed  beneath  the  ray  path  plot. 

D.   SUBROUTINE  RUNK 

Subroutine  RUNK  utilizes  a  fourth-order  Runge-Kutta 
numerical  method  to  integrate  the  internal  wave  equation  over 
the  water  column.   The  endpoint  of  the  integration  is  refined 
by  a  Newton— Raph son  iteration  method  until  the  endpoint  is 
l/200  of  the  maximum  value  of  the  eigenf unction ,  V,  determined 
during  the  integration  scheme.   For  each  selected  frequency, 
the  first  approximation  for  the  initial  eigenvalue  is 

K,  =7=  (45) 

\/gh 

The  approximation  is  refined  until  the  wave  number  correspond- 
ing to  the  first  mode  is  determined.   The  approximation  for 
the  second  mode  is  twice  the  value  of  K-^ .   The  initial  esti- 
mate of  the  remaining  eight  eigenvalues  is  the  previous 
eigenvalue  plus  the  difference  between  the  previous  two 
eigenvalues . 

After  the  first  three  eigenvalues  are  determined,  the 
program  plots  their  modal  structure.  The  remaining  seven 
wave  numbers  are  then  listed. 
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As  the  selected  frequencies  within  the  inertio-gravity 
band  increase  in  magnitude,  the  boundary  constraints,  though 
analytically  well-defined,  weaken  exponentially.   The  higher 
frequencies  are  greater  than  the  incremental  Vaisala  frequency 
for  a  greater  percentage  of  column  depth  and  it  becomes 
increasingly  difficult  for  the  numerical  program  to  find 
eigenvalues.   RUNK  terminates  its  search  for  the  next  eigen- 
values by  (l)  stopping  after  20  refinements  by  the  Newton- 

Raphson  iteration  process,  (2)  restricting  the  magnitude  of 

?  0 
the  eigenf unction  to  10    cm/sec,  and  (3)  stopping  if  the 

solution  converges  to  a  mode  other  than  the  one  searched  for. 

Also,  as  the  frequency  within  the  band  increases,  V 

becomes  more  sensitive  to  changes  in  k  when  seeking  higher 

modes.   For  situations  where  a  minute  change  in  k  causes 

abrupt  responses  in  end  values  of  V,  the  eigenvalue  may  be 

only  refined  to  a  certain  degree.   These  eigenvalues  are 

indicated  as  approximations. 
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IV.   TESTING  THE  PROGRAM 

A.   LINEAR  INPUT 

To  check  numerical  formulation,  coding,  and  accuracy,  two 
test  cases  were  used  whose  analytic  solutions  were  possible. 
The  first  test  was  a  linear  variation  of  depth  from  1.023 
grams/cc  at  the  surface  to  1.033  grams/cc  at  the  bottom  depth 
of  1000  meters.   The  Vaisala  frequency  for  this  input  was 
constant  and  equal  to: 

N 


r 


h^rrr     X  KT3  =9.77X  lO^sec"1.  (46) 

1.025    614/cc 

The  plot  by  Lhe  program  of  the  density  matches  that  of  the 
linear  variation.   The  Vaisala  frequency  determined  numeri- 
cally matches  that  determined  analytically.   Table  I  compares 
the  analytical  and  numerical  results  of  the  linear  variation 
of  density  for  frequency  0.00684  cycles/min.   Figure  4 
compares  the  first  three  modal  structures.   Analytically,  the 
normal  modes  are 


2   ^2  1/2 

m   n< 

The  eigenf unction  is: 


Km=(^f)    (f).  (47) 

N  -  f 


p  (£§)  sin  (fz)  (48) 


ex 

g 


B.   THREE-LAYER  INPUT 

The  second  input  was  a  density  strata  that  was  constant 
at  1.023  grams/cc  from  the  surface  to  300  meters,  varied 
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Table    I.  A   comparison    of    the    first    ten   normal    modes 
determined    for    frequency   0.00684    cycles/min: 

Analytical  Numerical  jo   Error 

k,  3.613 x 10~5  3.629 x  10~5  0.45 

k2  7.225 x 10~5  7.259 x 10" 5  0.47 

k3  1.083  xlO"4  1.088  xlO"4  0.27 

k4  1.445  xlO"4  1.451  xlO"4  0.42 

k5  1.806  xlO"4  1.814  xlO"4  0.44 

k.  2.167  xlO-4  2.177  x  10"4  0.46 

6  -4  -4 

k7  2.529x10^  2.540x10  0.43 

kg  2.890 x 10~4  2.903 x 10"4  0.45 

k„  3.251 x 10~4  3.266 x 10"4  0.46 


'9 


,-4  .    ,_       ,  „-4 


k1Q  3.613x10  3.629x10  0.44 


24 


-3 


VERTICAL      VELOCITY     (Cm/sEC     X    10       ) 


-  20 

-10 

0 

10 

20 

30 

I 

1 

i 

I 

~r 

I 

—  200 

i 
i 

—  400 

DC 


I 


0l 
UJ 

a 


—  60  0 


■  1000 


A— ANALYTICAL     (^N  O  RMA  L  IZE  D  ) 
X—  NUMERICAL  ' 


FIGURE  4.   Comparison  of  Eigenf unctions  for  First 
3  Modes  of  Frequency  0.00684 
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linearly  from  300  meters  to  a  value  of  1.027  grams/cc  at  700 
meters  ,  and  then  remained  invariant  from  700  meters  to  the 
bottom  at  a  depth  of  1000  meters.   The  corresponding  value 
of  the  Vaisala  frequency  was  then  zero  in  the  upper  and  lower 
isopycnal  layers  and  equal  to  9.8x10    in  the  middle  region. 

The  density  plot  by  the  program  matched  the  analytic 
solution  to  the  third  decimal  place.   Figure  5  illustrates 
the  comparison  between  the  analytic  solution  and  the  computer 
solution  for  the  Vaisala  frequency,  and  illustrates  how  the 
interpolation  scheme  handles  discontinuities  in  the  density 
slope . 
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FIGURE  5.   Comparison  of  Vaisala  Frequency  Plots 
by  Analytical  and  Numerical  Methods 
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V.   APPLICATION  OF  THE  PROGRAM  TO  OCEANOGRAPHIC  DATA 

The  program  utilizes  density  strata  from  Nansen  casts. 
It  can  process  several  stations.   IDATA,  the  number  of  sta- 
tions to  be  studied,  is  the  first  card  to  be  typed.   It 
follows  the  alphanumer ics  that  are  read  in  and  serves  in  the 
plotting  portions  of  the  program.   IDATA  is  read  in  as  an 
integer,  right-oriented  to  column  ten  of  the  IBM  card.   The 
second  card  to  be  typed  up  consists  of  three  figures.   The 
first  is  M,  the  number  of  bottles  at  the  station,  read  in  as 
an  integer,  right-oriented  to  column  10.   The  second  is  ZMAX, 
the  depth  at  which  the  last  bottle  was  taken,  read  in  the 
format  F10.1,  r  ight-or  ien  led  Lo  column  20.   The  third  figure 
is  RLAT ,  the  latitude  of  the  station,  that  is  read  in  using 
format  F10.1  right-oriented  to  column  30. 

The  next  cards  to  be  typed  are  those  representing  the 
density  strata.   These  are  read  in  as  four  pairs  of  numbers 
for  each  card.   Each  number  has  a  field  width  of  ten  and  is 
right-oriented.   The  first  number  in  each  pair  represents 
station  depth  and  has  a  format  of  F10.1.   The  second  of  the 
pair  has  a  format  F10.5  and  represents  the  station  density. 
These  four  pairs  of  numbers,  representing  80  columns  per  card, 
are  entered  until  all  the  station  data  has  been  read  in. 

The  above  represent  all  the  cards  needed  to  be  typed  for 
each  station.   From  experience  gained  using  the  IBM  360 
computer  to  process  ocean  stations  taken  in  mid-August,  1970, 
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off  Monterey,  California,  it  appears  that  two  minutes  are 
required  to  process  shallower  stations  (around  200  meters) 
and  four  minutes  to  process  deeper  stations  (around  2000 
meters )  . 

The  actual  printed  output  for  each  station  appears  as 
follows.   On  the  first  page  of  printed  output  appears: 
(l)  the  number  of  station  depths,  (2)  the  lowest  bottle  depth, 
(3)  a  printout  of  the  station  depths  and  densities  in  columnar 
form,  (4)  the  inertial  frequency  computed  for  the  station, 
(5)  the  maximum  Vaisala  frequency  computed  for  the  station, 
and  (6)  the  ten  frequencies  within  the  inertio-gravity  band 
selected  for  consideration. 

The  second  page  of  printed  output  plots  the  Vaisala 
frequency  and  density  versus  depth. 

The  next  ten  pages  plot  the  ray  path  for  each  frequency 
versus  depth.   At  the  bottom  of  each  ray  plot,  the  amount  of 
time  that  the  frequency  is  greater  than  the  Vaisala  frequency 
is  indicated. 

The  remaining  twenty  pages  of  output  display  information 
on  normal  modes.   Corresponding  to  each  of  the  ten  selected 
frequencies  are  two  pages.   The  first  page  lists  the  first 
three  normal  modes  and  plots  these  modes  versus  depth.   The 
second  page  lists  the  remaining  seven  modes. 

Two  stations,  one  shallow  and  one  deep,  were  chosen  to 
illustrate  the  output  from  the  program.   The  results  for  both 
stations  are  summarized  by  the  following  tables  and  figures. 
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Table  II  lists  the  information  for  the  deep  oceanographic 
station  displayed  on  the  first  page  of  computer  output.   This 
station  is  representative  of  an  oceanic  internal  wave  regime. 
Figure  6  shows  the  plot  of  the  density  variation  and  the 
Vaisala  frequency  for  the  deep  station.   The  Vaisala  plot 
follows  the  pattern  of  the  density  profile. 

Figure  7  depicts  the  ray  paths  for  the  first  three  select- 
ed frequencies.   The  first  frequency  selected,  0.01371 
cycles/min,  never  exceeds  the  incremental  Vaisala  frequency. 
Thus,  theoretically,  its  solution  is  always  oscillatory  and 
it  travels  horizontally  nearly  30  km.   The  second  frequency 
selected,  0.12969  cycles/min,  exceeds  the  Vaisala  frequency 
in  the  region  of  the  water  column  below  a  depth  of  970  meters. 
Thus  it  exponentially  'decays'  in  this  region  and  its  ray 
path  over  this  portion  of  the  water  column  is  vertical. 
Figures  7  and  8  illustrate  that,  with  increasing  frequency, 
the  horizontal  distance  covered  by  the  ray  path  decreases 
and  a  greater  portion  is  in  the  exponential  region.   The  ray 
paths  for  the  last  four  selected  frequencies  were  not  shown 
since  their  ray  is  nearly  vertical  from  top  to  bottom. 

Figure  9  plots  the  first  three  modes  for  frequency  0.01371 
cycles/min.   This  was  the  only  frequency  for  which  the  first 
three  mornal  modes  could  be  positively  identified.   Besides 
gaining  a  better  understanding  of  the  actual  appearance  of 
the  modal  structures,  Figure  9  suggests  optimum  placement  of 
a  thermocline  follower  for  mode  detection.   For  example,  the 
follower  would  best  detect  mode  three  at  a  depth  of  about 
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TABLE  II.   Data  dis 

play 

for  deep  station. 

Incremental  de 

pth 

Incremental  density 

0.0 

1.02498 

10.0 

1.02505 

20.0 

1.02522 

30.0 

1.02565 

50.0 

1.02588 

75.0 

1.02611 

100.0 

1.02624 

150.0 

1.02646 

200.0 

1.02656 

300.0 

1.02679 

400.0 

1.02681 

500.0 

1.02702 

600.0 

1.02711 

800.0 

1.02730 

1000.0 

1.02743 

1200.0 

1.02750 

1500.0 

1.02759 

2000.0 

1.02770 

2500.0 

i ;0?774 

3000.0 

1.02776 

Frequency  numt 

>er 

Cycles/min 

inertial 

0.00082 

1 

0.01371 

2 

0.12969 

3 

0.25855 

4 

0.38741 

5 

0.51628 

6 

0.64514 

7 

0.77401 

8 

0.90287 

9 

1.03173 

10 

1.16059 

maximum  Vaisala 

1.28946 
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2000  meters.   Table  III  lists  the  various  modes  determined 
for  the  selected  frequencies  for  the  deep  oceanographic 
station . 

Table  IV  displays  the  first  page  of  printed  output  for 
the  shallow  station  representative  of  the  regime  in  a  large 
bay.   Figure  10  illustrates  the  density  variation  and  the 
Vaisala  frequency  as  a  function  of  depth.   As  in  the  deep 
station,  the  first  selected  frequency  of  the  shallow  station 
has  an  oscillatory  solution  throughout  most  of  the  water 
column.   Its  ray  path  travels  horizontally  nearly  2500  meters. 
Its  exponential  region  may  be  seen  as  a  vertical  path  at  a 
distance  of  1955  meters.   Figure  11  shows  the  effect  of  an 
increase  in  frequency  without  a  considerable  increase  in  the 
effect  of  exponential  decay.   The  second  frequency  selected, 
0.15818  cycles/min.,  has  an  exponential  solution  over  the 
water  column  in  nearly  the  same  region  as  the  first.   However, 
frequency  2  only  travels  a  horizontal  distance  of  230  meters. 
Figure  12  depicts  the  decrease  in  ray  path  as  the  frequency 
increases . 

Figure  13  depicts  the  first  three  normal  modes  for 
frequency  0.01656  cycles/min.   Figure  14  illustrates  the 
first  three  normal  modes  for  frequency  0.31554  cycles/min. 
Figure  14  illustrates  well  the  manner  in  which  the  modal 
structure  deteriorates  with  increased  frequency.   Table  V 
lists  the  modes  determined  for  selected  frequencies  at  the 
shallow  station. 
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Table  III.   Normal  modes  for  selected  frequencies  at 
the  deep  station. 

0.01371  cycles/min.         0.12969  cycles/min 


1. 

0.00011736 

1. 

0.0018469 

2. 

0.00023344 

2. 

0.0034051 

3. 

0.00032185 

4. 

0.00043948 

5. 

0.00055503 

6. 

0.00063828 

7. 

0.00079211 

8. 

0.00085791 

9. 

0.00096286 

10. 

0.00011641 
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Table  IV.   Data  display  for  shallow  station. 
Incremental  depth     Incremental  density 


0.0 

1.02535 

5.0 

1.02537 

10.0 

1.02543 

15.0 

1.02551 

20.0 

1.02570 

25.0 

1.02580 

30.0 

1.02606 

40.0 

1.02610 

50.0 

1.02614 

60.0 

1.02614 

75.0 

1.02616 

100.0 

1.02619 

Frequency  number 

Cycles/min 

inertial 

0.00083 

1 

0.01656 

2 

0.15818 

3 

0.31554 

4 

0.47289 

5 

0.63024 

6 

0.78760 

7 

0.94495 

8 

1.10231 

9 

1.25966 

10 

1.41702 

maximum  Vaisala 

1.57437 
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Table  V.   Normal  modes  for  selected  frequencies  at  the 
shallow  station. 


0.01656  cycles/min 


1. 

0.00087701 

2. 

0.0025324 

3. 

0.0035674 

4. 

0.0050806 

5. 

0.0055217 

o.: 

31554  cycles/min. 

1. 

0.018642 

2. 

0.056013 

3. 

0.10871 

4. 

0.14614 

5. 

0.17507 

6. 

0.20301 

7. 

0.23839 

8. 

0.28283 

0.15818  cycles/min 


10.   0.37399 

0.47289  cycles/min 


1. 

0.0086086 

2. 

0.025461 

3. 

0.041537 

4. 

0.051653 

5. 

0.069448 

6. 

0.074546 

7. 

0.092209 

8. 

0.10860 

9. 

0.12385 

10. 

0.13605 

0.63024    cvcles/min 


1. 

0.032306 

1.   0.053539 

2. 

0.094660 

2.   0.14805 

3. 

0.19395 

4. 

0.26040 

0.78760  cycles/min. 

5. 

0.33316 

6. 

0.40811 

1.   0.092463 

7. 

0.48565 

2.   0.23132 

8. 

0.55884 

9. 

0.63999 

.0. 

0.72820 
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VI.   CONCLUSIONS 

The  numerical  program  developed  in  this  thesis  has  derived 
a  technique  that  links  the  theoretical  application  of  internal 
waves  to  real  data.   Real  solutions  are  determined  for  actual 
density  profiles.   The  program  can  yield  a  useful  description 
of  the  internal  wave  regime  within  a  given  area.   From 
representative  samples  of  stations  investigated,  four  points 
are  apparent: 

(l),  that  with  greater  depth  the  effect  of  the  boundary 
conditions  weakens; 

(2),  that  for  low  frequencies  that  exceed  the  Vaisala 
frequency  over  only  a  small  percentage  of  the  water  column, 
the  ray  path  approaches  the  horizontal,  while  for  higher 
frequencies  the  ray  path  approaches  the  vertical; 
(3),  that  the  plot  of  the  lowest  three  modes  demonstrates 
the  manner  in  which,  with  increasing  frequency,  the  effect 
of  the  boundary  conditions  weakens  as  the  oscillatory 
regions  give  way  to  exponential  "decay"  regions; 
(4),  that  for  the  stations  investigated  by  this  program, 
only  those  frequencies  within  the  lower  quarter  of  the 
inertio-gravity  frequency  band  could  be  solved  for  normal 
modes  with  any  degree  of  success. 
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